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Abstract 

In  this  paper  we  discuss  a  general  methodology  for  estimating  the  dis¬ 
tribution  of  individual  growth  rates  in  a  size-structured  population  using 
aggregate  population  data.  The  method,  for  which  rigorous  theoretical 
formulations  have  been  developed,  is  presented  in  the  context  of  an  in¬ 
verse  problem  methodology  and  its  use  is  illustrated  with  application  to 
mosquitofish,  Gambusia  affinis,  population  in  rice  fields. 


1  Introduction 

In  this  paper,  we  present  results  using  inverse  problem  techniques  for  estimation 
of  growth  distribution  in  size-structured  population  models  using  aggregate  pop¬ 
ulation  data.  The  models  employed  here  are  based  on  ideas  initially  discussed  in 
[BBKW],  which  entail  models  wherein  growth  rates  may  vary  across  individuals 
of  the  population  as  well  as  with  size  and  time. 

These  models  are  in  contrast  to  the  usual  stochastic  partial  differential  equa¬ 
tion  models  as  described  for  example  in  [FL1],[FL2],[FS]  and  in  the  next  section. 
Although  they  are  not  stochastic  in  the  usual  sense,  they  are  probabilistic  in 
that  one  models  growth,  mortality,  etc.  via  probability  distributions  across  the 
population. 

Simulation  studies  were  presented  in  [BBKW]  to  demonstrate  that  such  ideas 
could  lead  to  population  densities  that  exhibit  dispersion  and  bimodality.  Rig¬ 
orous  theoretical  developments  of  the  associated  inverse  problem  technique  and 
initial  illustrations  with  computational  examples  were  given  in  [BF],  [FI]  and 
[F2] .  A  survey  of  results  and  other  references  can  be  found  in  [B],  The  general 
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philosophy  underlying  our  approach  of  using  aggregate  population  data  to  esti¬ 
mate  individual  rates  is  given  (along  with  an  application  to  susceptibility  and 
vaccination  efficiency  in  populations)  in  [BFZ], 

As  detailed  in  earlier  references,  our  efforts  on  such  problems  were  initiated 
in  collaboration  with  marine  biologists  (specifically  Lou  Botsford  and  his  col¬ 
leagues)  in  studies  related  to  the  introduction  of  mosquitofish  into  rice  fields,  in 
attempts  to  control  mosquito  populations  without  chemicals.  To  establish  viable 
control  strategies,  it  is  very  important  to  have  quantitative  models  which  pre¬ 
dict  accurately  the  evolution  of  the  populations.  In  the  paper  [BVWLKRC],  the 
mosquitofish  population  was  modeled  using  the  Sinko-Streifer  (also  called  the 
McKendrick-Von  Foerster)  model  for  size-structured  population  density  evolu¬ 
tion.  As  we  shall  discuss  below,  this  leads  to  a  number  of  conceptual  difficulties. 

2  Modeling  Philosophy 

There  is  a  huge  literature  (we  won’t  even  attempt  to  cite  these  here)  on  the  mod¬ 
eling  of  populations  in  the  quantitatively  oriented  biological  research  literature. 
The  efforts  range  from  modeling  population  growth,  death,  mortality  in  insects, 
marine  life,  plants  and  animals  to  gene  frequency  mutation  and  drift,  as  well  as 
treatments  of  susceptibility  to  disease  in  vaccinated  populations  of  humans.  In 
such  a  diverse  scientific  literature  it  is  often  difficult  to  discern  an  underlying 
commonality.  However,  there  are  some  usual  features  in  much  of  the  literature. 
The  modeling  often  relies  on  data  about  certain  observables,  while  it  is  knowl¬ 
edge  about  other  nonobservable  parameters  that  is  of  interest  to  investigators. 
For  example,  one  usually  has  total  population  (aggregate  densities)  counts  in 
models  where  growth,  mortality  and  migration  rates  of  a  typical  individual  are 
of  most  interest.  Or  one  may  have  aggregate  data  on  percent  of  a  population 
vaccinated,  and  numbers  of  those  who  fall  ill  and  recover  (or  not)  in  studies  of 
disease  prophylaxis.  In  the  case  of  genetic  studies  one  has  data  on  phenotypes, 
where  it  is  gene  frequencies  and  their  changes  that  one  wishes  to  understand. 

A  second  feature  of  most  modeling  attempts  is  the  presence  of  uncertainty. 
This  may  arise  in  the  model  itself,  in  parameters  in  the  model,  in  distribution 
over  the  population  of  unobservable  traits  or  characteristics.  A  conceptually  im¬ 
portant  question  is  how  to  treat  these  uncertainties  when  attempting  to  extract 
the  maximum  information  about  the  population  from  the  data.  The  standard 
approach  to  treat  probabilistic  or  stochastic  aspects  of  population  is  through 
the  use  of  stochastic  differential  equations.  A  brief  review  of  examples  early  on 
in  this  approach  was  given  by  Fleming  in  [FL1],  where  he  discussed  distributed 
parameter  or  partial  differential  equation  models  for  geographically-structured 
(including  size  or  age-structured)  populations.  He  also  discussed  population 
genetics  models  (equations  for  the  mean  and  the  covariance  of  gene  frequency 
in  spatially  migrating  populations  —  see  [FL2],[FS]).  Our  focus  here  will  be 
on  size-structured  population  models  and  our  use  of  rate  distribution  models  as 
an  alternative  to  stochastic  differential  equations.  To  put  this  in  context  of  the 
more  standard  approach,  we  introduce  and  discuss  briefly  the  Fokker-Planck 
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size-structured  model. 

The  Fokker-Planck  equation  is  the  basis  of  a  stochastic  size-structured  model 
[B1],[BTW]  which  is  based  on  the  assumption  that  movement  from  one  size 
class  to  another  can  be  described  by  a  Markov  process.  The  “physiological  age” 
interpretation  of  the  Fokker-Planck  equation  was  first  suggested  by  Weiss  [W];  a 
careful  derivation  based  on  the  paradigm  of  Brownian  motion  of  particles  which 
is  applicable  to  growth  processes  is  given  in  [O].  The  Fokker-Planck  equation, 
under  the  assumption  of  a  Markov  transition  for  the  growth  process,  is 

a--)  +  x)u(t ,  a;))  =  x)u(t,  x)),  (1) 

where  n(t.  x)  is  the  population  density  at  time  t  and  size  x  and  the  moments 
are  given  by 

Mj(t,x)=  lim  —  /  (y  -  x)Jp(t,x]t  +  At,y)dy.  (2) 

At^O  /At  J_OD 


The  function  p(t,x]t  +  At,y)  is  the  probability  density  for  the  transition 
from  size  x  at  time  t  to  size  y  at  time  t  +  At]  i.e. ,  p(t,  x;  t  +  At,y)Ax  is  the 
probability  that  members  in  the  size  interval  [x,x  +  Aa]  at  time  t  will  move  to 
size  y  at  time  t  +  At.  The  moments  M\,  M2  have  probabilistic  interpretations: 
M  i  is  the  mean  (or  first  moment)  of  the  rate  of  increase  in  size 


M.H: 


=  lim  £  ■ 

At^O 


s(t  +  At)  —  a :(t) 
At 


}• 


where  ^-[A']  denotes  the  expected  value  of  a  random  variable  A',  while  M2  is 
the  second  moment  of  the  rate  of  increase  in  size 


M2{t,x) 


lim  £  { 

At^O  ( 


(x(t,  +  At)  —  x(t ))2  | 

At  J  ' 


Appropriate  boundary  conditions  must  be  formulated  for  (1).  Since 
M  |  u  —  ^  represents  the  population  flux,  we  have 


.Mi (t,  x)u(t,  x) 


—  (M2(t,x)u(t,x)) 


Mi(t,  x)u(t,  x) 


('Em 

/  K(t,  x)u(ta  x)da 
J  Xn 


—(M2(t,x)u(t,;x)) 


0 


(3) 

(4) 


while  the  initial  conditions  are  given  by 


u(t0,  x)  —  <F(a:).  (5) 

Here  xm  is  maximum  size,  x0  is  the  minimum  size,  and  K  is  the  fecundity 
rate  as  explained  below  in  our  discussion  of  the  Sinko-Streifer  model. 
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The  system  (1),  (3)-(5)  comprises  an  initial-boundary  value  problem  for 
a  transport-dominated  diffusion  process  (the  M  \  term  is  typically  dominant 
over  the  M2  term)  that  offers  significant  computational  challenges.  First,  the 
moments  Mi,  M2  almost  always  (for  the  applications  to  populations)  depend 
on  both  t  and  x  and  must  be  estimated;  the  expressions  (2)  are  not  helpful 
since,  in  general,  p  is  unknown.  Moreover,  even  simulation  with  (1),  (3)-(5) 
is  nontrivial.  Traditional  finite  difference  and  finite  element  methods  produce 
erroneous  oscillatory  solutions.  Thus  the  spline-based,  fixed  node  methods  first 
proposed  in  [Bl]  are  of  very  limited  interest.  However,  high  promise  can  be  found 
[BTW]  in  a  moving  node  finite  element  technique  first  suggested  in  [H],  Even 
so,  estimation  or  inverse  problems  built  on  Fokker-Planck  models  are  extremely 
different  computationally  and  have  enjoyed  limited  use  in  the  literature. 

Our  efforts  to  treat  both  unobservable  individual  parameters  as  well  as  a 
certain  amount  of  stochasticity  is  built  upon  the  Sinko-Streifer  model  and  thus 
we  briefly  outline  its  features. 

The  classical  Sinko-Streifer  model,  henceforth  referred  to  simply  as  (SS),  is 
given  by 

x)  +  x)v(t,  x))  =  -p(t,  x)v(t,  .!•)  t  >  t0  , 

t;(0,  *)  =  $(*)  xo  <  x  <  xm  , 

(6) 

g(t,x0)v(t,x 0)  =  K(t,Z)v(t,Z)d£ 

g{t,xm )  =  0 


and  simulates  the  time  evolution  of  a  population  with  respect  to  the  size  x  of 
the  individuals.  In  our  use  of  this  model,  the  parameter  x  £  [xo,xm]  denotes 
the  length  of  the  fish,  and  the  function  v  represents  the  size  density  function  so 
that 

N(t)  =  I  v(t,  x)dx  (7) 

J  a 

is  the  number  of  fish  in  the  population  at  time  t  whose  size  is  between  a  and  b. 
The  function  g  represents  the  growth  rate  of  the  individuals, 


dx(i) 

dt 


g(t,x)  , 


so  that  in  this  simple  model  all  individuals  have  the  same  growth  rate.  The  last 
condition  in  (6),  g(t,xm )  =  0,  indicates  that  growth  ceases  when  individuals 
reach  maximum  size  xm.  The  initial  condition  $  describes  the  initial  size  density 
for  the  population.  The  function  K  is  a  fecundity  kernel,  and  is  used  to  express 
the  recruitment  rate  R(t,v )  =  fx™  K(t,!;)v(t,!;)dt;.  The  mortality  rate  is  given 
by  /i. 

Difficulties  arise,  however,  in  applying  the  model  (SS)  to  observed  data  since 
the  data  often  exhibit  features  that  the  model  cannot  simulate  or  predict.  One 
such  phenomenon  is  a  dispersion  in  size  as  time  progresses.  Another  is  that  the 
population  begins  with  a  unimodal  density  and  in  time  develops  into  one  with 
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a  bimodal  density.  Both  of  these  features  are  present  in  the  mosquitofish  data 
(see  [BBKW],  [BVWLKRC],  and  the  data  plot  in  Figure  1). 

This  dispersion  in  size  and  the  transition  from  a  unimodal  to  a  bimodal 
density  is  qualitatively  inconsistent  with  solutions  of  the  Sinko-Streifer  model 
under  biologically  feasible  assumptions  on  the  model.  We  are  able  to  capture 
both  phenomena  with  growth  rate  distribution  (GRD)  models. 

2.1  GRD  Models  with  Identical  Initial  Size  Densities 

Here  we  use  a  modification  of  the  model  (SS)  (in  actuality,  a  continuum  of  (SS) 
across  families  of  growth  rates)  which  exhibits  features  present  in  the  data.  We 
use  as  the  aggregate  population  density  (APD)  the  function 

u(t,x)=  /  v(t,  x;  g)dP(g);  (8) 

JG 

where  G  is  a  collection  of  growth  rates  and  P  is  a  probability  measure  on  G. 
This  approach  was  first  suggested  in  [BBKW].  As  detailed  in  [BF],  the  resulting 
model  is  rich  enough  to  exhibit  the  phenomena  of  interest,  namely,  dispersion 
and  development  of  two  modes  from  one.  As  also  explained  in  [BF],  measure 
theoretic  results  allow  us  to  approximate  the  continuum  measure  in  (8)  by  a 
discrete  measure  corresponding  to  a  finite  dimensional  set  of  growth  rates. 

Briefly,  we  assume  that  we  have  a  family  G  =  {gjk},  j  =  1,  •  •  ■  ,M\,  k  = 
!,■■■,  M2  of  individual  growth  rates;  individuals  are  grouped  into  the  same 
subpopulation  if  they  have  the  same  growth  rate.  The  subpopulations  grow 
according  to  the  (SS)  with  g  =  gjk  in  (6).  Here  we  take  gjk{%)  =  rj( 7*,  —  j:). 
H  =  0  and  K  =  0.  Note,  however,  that  a  generalization  to  allow  distribution  of 
mortality  and  fecundity  over  subpopulations  is  readily  achieved. 

The  size  density  for  the  subpopulation  jk  is  given  by  v(t,x;  gjk),  and  the 
aggregate  population  density  is  given  by 

u(t,  x)  =  ^  x\  9jk)Pjk ,  (9) 

j,k 

where  pjk  is  the  probability  of  an  individual  being  in  the  jk  subpopulation. 
Note  that  this  formulation  embodies  the  assumption  that  each  subpopulation 
has  the  same  initial  size  density,  i.e. ,  u(0 ,  x',gjk)  =  $(*)  for  all  j,k. 

We  used  the  above  outlined  formulation  in  a  least  squares  inverse  problem 
involving  fitting  several  sets  of  field  data  {u(t,x)}.  The  probabilities  {pjk}  for 
1  <  j  <  M 1  and  1  <  k  <  M->  in  (9)  were  estimated  by  solving  the  following 
inverse  problem: 
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min 

P  ePM(G) 


j(  p) 


N  r1  •> 

'y  '  /  i)  :  •’  ■  flj k  )pj h.(f  n  ,  X )]  (10) 

n  =  l  x°  j,k 
N  rx  i 

E  «E  ^  (^n  3  ^5  9j  k  ^Pj  k  ] 

n  =  1^a7°  j,ic 

2«( /„  ■■)-)  y:  y{t„,x-,  gjk)pjk  +  [u(tn,x)]2d,x. 
j,k 


Here  PM(G)  is  a  finite  dimensional  approximation  to  the  probability  measure 
space  P(G )  defined  by 

pm(g )  =  {Pe  P(G)\p  =  Y'Pjktg^,  £>,-*  =  1} 

j,k  j,k 

where  M  =  M\  x  M2  and  5gjk  is  the  Dirac  measure  with  an  atom  at  gjk.  We 
denote  by  p  the  array  that  contains  pjk,  1  <  j  <  Mi,  1  <  k  <  M2,  and  set 

N  t*  1 

G’jklm  —  E  v(tn  ,  X]  yj k)  t.( I, f  ,  X]  (Jhn  )dx 

n  =  1  Jx« 


N  ^x1 

bjk  =  -  E  /  u(tn,x)v(tn,x-,  cjjk)dx 

n  =  1  Jx« 

N  r%  1 

c  =  /  [u(tn ,  x)]2c?a:  . 

n  =  l  •'*<» 

Moreover,  we  define  A  and  b  to  be  the  arrays  that  contain  a.jkim  and  bjk 
respectively,  for  1  <  j.  f  <  M\  and  1  <  k,m  <  M2- 

The  inverse  problem  (10)  now  becomes  a  quadratic  programming  problem  in 
which  one  minimizes  pTAp  +  2pTb  +  c  over  PM(G).  To  include  the  equality 
constraint,  ^ pjk  =  1,  we  introduce  the  Lagrange  multiplier  A  and  solve  the 
unconstrained  problem  minF(p,A),  where 

F( P,  A)  =  pTAp  +  2pTb  +  c  +  K'y^Pjk  ~  1]  •  (11) 

j,k 


We  used  the  “method  of  characteristics”  based  techniques  discussed  in  [IKP] 
to  solve  (6)  for  the  densities  v(t,x;  gjk),  which  were  calculated  in  parallel  using 
PVM  software  on  six  IBM  RISC  6000  workstations.  These  densities  were  then 
used  with  the  IMSL  subroutine  DQP  ROG  to  solve  for  the  optimal  solution 
(p*,A*).  Finally,  the  aggregate  population  density  u(t,x )  was  determined  by 
substituting  p*k  into  (9). 
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In  the  above  parameter  estimation  problem,  for  each  subgroup  jk,  a  propor¬ 
tioned  initial  size  density  <&(x)pjk  was  assumed.  The  Day  195  data  (see  Figure 
1)  was  chosen  so  that 

^v(t0,x;  gjk)pjk  =  ^2  ®(x)Pjk ■ 

j,k  j,k 

Note  that  since  the  first  data  set  is  used  here  as  the  initial  size  density,  it  cannot 
be  used  in  solving  the  inverse  problem.  This  leaves  one  less  data  set  available 
for  use  in  estimating  the  optimal  parameter  set. 

Example  1:  We  present  in  the  following  an  example  which  is  typical  of  the 
field-data  fitting  results.  The  fish  were  stocked  on  June  28,  1982  into  four  rice 
paddies  with  parallel  water  flow.  Each  paddy  had  an  outflow  trap  to  measure 
emigration.  Measurements  were  taken  weekly  (two  paddies  one  day  and  two  the 
next)  with  the  use  of  fifteen  traps  per  paddy.  The  total  number  of  fish  caught 
was  greater  than  or  equal  to  the  number  actually  measured. 

The  size  distribution  frequency  for  size  class  i  is  defined  as  /,;  =  nmti/Nm , 
where  nm,i  is  the  number  of  fish  measured  in  size  class  i,  and  Nm  is  the  total 
number  of  fish  measured. 

The  total  population  was  divided  into  five  hundred  and  twelve  subpopula¬ 
tions  with 

TV  =  0.2+  1/31  *4.8*  (j  -  1)  j  =  1,2,...,  32 
7i  =  16/38,  72  =  22/38,  73  =  24/38, 

7i t  =  16/38  +  1/15  *22/38*  {k  -  1)  k  =  4, ...,  16. 

The  range  for  rj  was  chosen  arbitrarily,  while  the  7*,  were  based  on  the  field 
data.  The  choice  for  71  =  16/38  was  made  by  inspecting  the  data  from  Day 
195.  At  x  =  16mm,  a  significant  decrease  in  the  number  of  fish  occurs,  which  we 
interpreted  as  the  smallest  possible  maximum  length  for  some  subpopulation. 
We  obtained  71  =  16/38  and  7i@  =  1.0  after  normalizing  against  the  largest 
possible  length,  xi  =  38 mm,  for  the  total  population.  The  initial  size  density 
$(*)  (for  all  subpopulations)  was  approximated  by  interpolating  the  Day  195 
data. 

The  computed  results  manifest  dispersion  and  bimodality  as  we  expected 
(see  Figure  1),  with  a  residual  of  «7i  =  3.7626  x  10-4.  Figures  2  and  3  illustrate 
the  probability  densities  and  distribution  respectively  as  a  function  of  r  and  7. 

Although  the  computed  results  do  provide  dispersion  and  bimodality,  the 
mismatches  between  the  computed  solutions  and  the  field  data  are  still  signifi¬ 
cant.  Furthermore,  the  discrepancies  did  not  diminish  as  fast  as  hoped  when  we 
increased  the  number  of  subpopulations,  M  =  M 1  x  Mn.  The  idea  of  varying 
the  initial  size  density  among  the  subgroups  became  appealing. 

2.2  GRD  Models  with  Subpopulation  Dependent  Initial 
Size  Densities 

A  further  refinement  (and  one  which  is  appropriate  for  the  mosquitofish  popu¬ 
lation  presented  here)  entails  a  parameterization  of  the  initial  density 
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<&(*)  =  with  4>jk{x)  =  J2iXe(x)atjik.  Here  Xt{x)  is  the  character¬ 

istic  function  of  the  interval  [x, i,xi+i)  corresponding  to  a  partition  0  =  x0  = 
Xi  <  xn  <  ■  ■  ■  <  xm3+i  =  1,  and  we  have  normalized  the  size  range  [x0,  xm]  to 
[0,1].  If  one  solves  (6)  with  g  =  gjk,*&  —  Xli  then  v(t,  x\ gjk,  Xt)  represents  the 
density  for  individuals  with  growth  rate  gjk  that  initially  have  size  structure 
density  given  by  xt- 

The  subpopulation  jk  density  is  given  by 

I'it-.xigjk)  =  '^2v(t,x;gjk,xi)atj,k 
t 

so  that  the  aggregate  population  density  is  then  given  by 

u(t,  x)  =  ^2  v(t,  X ;  gjk)pjk  =  X]  X]  v(t’  x’  9jk,Xi)Pjkai-jtk-  (12) 
i,k  t  j,k 


As  before,  we  require  the  constraint  JA  kPjk  —  1.  We  denote  ct  as  the  array 
containing  utj,k  for  1  <  j  <  Mi,  1  <  k  <  M2  and  1  <  £  <  M3. 

This  formulation  leads  to  a  new  inverse  least  squares  problem  with  objective 
function 


J(p,a) 


u(tn ,  x ) 


i(tn  ;  X )  | 


dx 


X ,  gjk ,  Xl) Pj  k  C'7' :  j ,  k 


(13) 


2 


u(t„,x) 


dx 


subject  to  k  pjk  =  1,  and  v(t ,  X]  gjk,  Xt)  being  the  solution  of  (6)  correspond¬ 
ing  to  g  =  gjk  and  $  —  Xt-  Minimizing  (13)  over  both  p  and  a  requires  solving 
a  nonlinear  programming  problem. 

An  alternative  and  simpler  approach  to  directly  solving  (13)  is  to  reduce  the 
nonlinear  programming  problem  back  to  a  quadratic  programming  problem. 
This  can  be  accomplished  by  eliminating  either  p  or  a  from  the  minimization 
process. 


2.2.1  Minimizing  Over  a  Only 

One  method  for  reducing  (13)  to  a  quadratic  programming  problem  is  to  remove 
p  from  the  parameter  estimation  problem  and  simply  minimize  over  a.  That 
is,  if  we  can  pre-determine  a  good  choice  for  p,  we  can  fix  the  probabilities  pjk 
in  (13)  and  reformulate  it  as  a  quadratic  programming  problem  with  objective 
function 

F(a)  =  aTAo  +  2aTb  +  c.  (14) 

Here  A  is  the  M \  M-> M3  x  M1M2M3  array  containing  the  elements 
N  ,1 

ajkt,qst  =  V  /  [v{t„,x-,gjk,xt)pjk\[v{tn,x-,gqs,xt)pqs]dx,  (15) 

n  =  0J° 
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b  is  the  vector  with  elements 


N  »i 

bjki  —  'y  '  I  j x}v(tn ,  x ,  gj k ,  Xj?)Pj kdx y  (16) 

n  =  0J° 


and 


N 

E 

n  =  0 


f 

Jo 


[w(tn,  a;)]2c?a;. 


(17) 


A  good  way  to  pre-determine  the  probabilities  pjk  is  to  first  solve  the  problem 
for  (11)  with  identical  initial  size  densities.  We  can  then  use  the  solutions  p*  as 
the  fixed  probabilities  in  the  arrays  A  and  b. 

A  major  advantage  of  this  approach  is  that  the  formulation  allows  for  differ¬ 
ent  initial  size  densities  across  subpopulations.  Moreover,  since  the  initial  size 
densities  are  not  estimated  directly  with  the  first  data  set,  all  of  the  data  sets 
may  be  used  in  solving  the  inverse  problem. 

As  before,  we  used  “method  of  characteristics”  based  techniques,  parallel 
computation,  and  the  IMSL  subroutine  DQP  ROG  to  solve  for  the  optimal  so¬ 
lution  o;* .  The  resulting  k  and  the  pre  determined  p*-k  were  then  substituted 
into  (12)  to  solve  for  u(t,x)  corresponding  to  (p*, a*). 

Example  2:  In  the  following  example,  we  use  the  1982  field  data  as  described 
in  Example  1.  First  the  probability  measures  pjk  were  obtained  by  optimiz¬ 
ing  (11),  where  the  total  population  was  divided  into  thirty-two  subpopulations 
according  to  the  following: 


rj  =  0.2  +  1/7  *  4.8  *{j-  1)  j  =  1, 2, . . . ,  8 
7i  =  16/38,  72  =  22/38,  73  =  24/38,  74  =  1. 


(18) 


As  in  Example  1,  the  Day  195  data  was  used  as  the  initial  size  distribution 
( .7: )  for  all  subgroups. 

The  computed  probabilities  pjk  were  then  used  in  solving  the  reformulated 
quadratic  programming  problem  for  (14)  with  subpopulation-dependent  initial 
size  densities.  The  entire  population  was  further  divided  into  a  total  of  six 
hundred  eight  subpopulations,  with  rj  and  7^  as  in  (18)  for  j  =  1,2,...,  8, 
k  =  1, 2,  3, 4,  and 


X,  =  [xt,xt+ 1)  =  [1/19  *(/  1),  1/19  *£)  £  =  1, . . . ,  19. 

The  intervals  each  correspond  to  a  2mm  size  class  as  seen  in  the  field  data. 

The  computed  weights  k  were  then  substituted  into  (12)  along  with 
the  precomputed  p*-k  to  obtain  the  aggregate  population  density  u(t,x )  (see 
Figure  4).  This  calculated  density  gave  a  better  fit  to  the  field  data  than  the  fit 
obtained  in  Example  1.  Figures  5  and  6  illustrate  the  probability  densities  and 
distribution  respectively  as  a  function  of  r  and  7. 

The  overall  residual  was  J2  =  3.08780x  10-4.  Since  the  parameter  estimation 
carried  out  in  Example  1  did  not  utilize  the  Day  195  data,  we  define  a  second 
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Length  of  Fish  (mm) 


Figure  4:  Computed  results  (-  .)  vs.  field  data  ( — )  for  Example  2. 


residual  J2  which  is  taken  only  over  Days  202,  209  and  216  of  the  data.  This 
provides  a  means  of  comparing  the  accuracy  of  Example  1  versus  Example  2. 

In  this  manner  we  obtained  J2  =  2.9305  x  1 0 —  5 ,  which  is  quite  an  improve¬ 
ment  over  J ■  =  3.7626  x  10-4. 

2.2.2  Iterative  Quadratic  Programming 

An  alternative  approach  for  simplifying  the  nonlinear  programming  problem 
for  (13)  is  an  iterative  process  involving  two  separate  quadratic  programming 
problems.  The  first  of  these  quadratic  programming  problems  is  formulated  by 
minimizing  (13)  over  pjk  while  ottj, k  is  held  constant,  and  the  second  minimizes 
over  while  pjk  is  held  constant.  These  quadratic  programming  problems 
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can  be  solved  in  succession,  creating  an  iterative  process  that  progressively 
improves  upon  p  and  a. 

Since  the  variables  otp-j^  are  associated  with  the  initial  size  densities  <j>jk,  we 
can  approximate  oipj  j,  using  the  initial  time  data.  That  is,  we  let 

Xl  —  ['/  ^ f’ +  I  )  ^  —  1)  •  •  •  )  44  3 , 

subdividing  the  size  classes  in  the  same  manner  as  the  data.  For  each  l ,  we  then 
set 

—  h  ( 0 ,  Xp  ) 

for  every  j,k. 

These  values  of  apj^  can  then  be  inserted  as  fixed  quantities  in  (13)  to  arrive 
at  a  quadratic  programming  problem  with  objective  function 

F(p,  X)  =  pTAp  +  2pTb  +  c  +  X[^2  Pjk  ~  1]-  (19) 

j,k 


Here  A  is  the  Mi  M2  x  Mi  M2  array  with  elements 
N  |  M:\  A/ 3 

djktm.  —  'y  '  I  ■  9j k:  Xi)o?/;j;/]  1  y  (jPrr:  ■.  Xi )  ; kni ] d-C , 

n=0^°  |=1  i= 1 

the  vector  b  contains  elements 

N  „i  M3 

bjk  -  y:  /  »(/,; ,  a--)[y^  u(tn  ,  x ;  gjk:  Xi)<*i;j  ,k]dx , 

n  =  0^°  i  =  l 

and 

A 

c  =  ai)]2ejsr, 

„=oJ» 

Again  we  used  the  “method  of  characteristics”  based  techniques,  paral¬ 
lel  computation  and  the  IMSL  subroutine  DQP  ROG  to  obtain  the  optimal 

(P'-A'). 

These  probability  measures  p*- k  may  now  in  turn  be  used  to  solve  a  quadratic 
programming  problem  over  ai.j.k-  In  fact,  by  holding  p  constant  in  the  original 
problem  for  (13),  we  have  reduced  the  nonlinear  programming  problem  to  the 
quadratic  programming  problem  involving  ( 14 )-( 17) .  From  this  we  obtain  the 
optimal  a* . 

This  process  can  be  repeated,  leading  to  the  following  iterative  method: 

1.  Partition  the  size  interval  into  subintervals  {xp,  xp+i},  £  —  1,2, ... ,  M3,  so 
that  the  subintervals  correspond  to  the  subdivisions  of  size  classes  in  the 
data.  For  each  £,  let  0$  ■  k  =  w(0,  xp). 

2.  (a)  For  /  =  0,1,...,  let  =  (i’p  j  j,  and  solve  the  quadratic  program¬ 

ming  problem  for  (19)  to  obtain  the  optimal  p^k. 
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(b)  With  pjk  —  pljk  for  all  j,  k,  solve  the  quadratic  programming  problem 
for  (14)  to  obtain  the  optimal  o/4)  ). . 

3.  Stop  iterating  when  |p!  +  1  — p*|  <  TOL\  and  |td+1(0,  x)  — «J(0,  x)\  <  TOLn, 
for  some  user-specified  tolerances  TOL\  and  TOL 2. 

Our  choice  of  the  above  stopping  criteria  was  motivated  by  the  structure 
of  our  model.  The  pjk  are  probabilities  associated  with  the  subpopulations 
gjk,  and  enter  the  model  through  (9).  We  would  like  |p®  —  p*+1|  to  be  within  a 
specified  tolerance.  An  additional  stopping  criterion  should  involve  the  variables 
ai;j,k,  which  are  weights  associated  with  the  parameterization  of  the  initial  size 
density  <h(a;)  =  •  k^2iXi(x)ai;j,k-  The  individual  themselves  need  not 

converge  upon  iteration,  but  we  would  like  the  initial  size  densities  to  converge. 
That  is,  we  would  like 

|ff  (0,::*)  —  if  +  (0.  .r)|  <  TOL2,  for  some  specified  tolerance  TOZ2. 

Example  3:  In  this  example  we  continue  to  use  the  1982  field  data  as  detailed 
in  Example  1.  To  initialize  the  iterative  process,  we  approximated  the  weights 
atj  k  by  the  initial  data  from  Day  195.  That  is,  we  let 

Xi  =  [xi,xt+ 1)  =  [1/19  *  ii  —  1),  1/19*1)  £=  1,...,19, 

and  for  each  £  we  set  a%j  k  =  w(0,  xi)  for  all  j,  k. 

We  divided  the  total  population  into  six  hundred  eight  subpopulations  with 

rj  =  0.3  +  1/3  *  4.7  *  (j  -  1)  j  =  1,  2, 3, 4 

71  =  16/38,  72  =  22/38,  73  =  24/38,  (20) 

7*  =  16/38+  1/7 *22/38 *{k  -  1)  jfe  =  4,..., 8. 

The  variables  0 rj  ■  k  were  then  inserted  into  (19)  and  the  quadratic  program¬ 
ming  problem  was  solved  for  the  probability  measures  p°-k.  The  corresponding 
aggregate  population  density  u°(t,x )  was  formed  by  substituting  •  k  and  p°-k 
into  (12). 

These  computed  solutions  u°(t,x )  did  not  fit  the  field  data  as  well  as  any 
of  the  previous  computed  solutions  (see  Figure  7),  returning  a  residual  of  J3  = 
0.0020,  and  a  residual  over  Days  202,  209  and  216  of  J3  =  0.0018. 

Continuing  the  iterative  process,  we  used  p°-k  to  optimize  (14)  for  a\.- k. 
This  led  to  a  new  set  of  computed  solutions  w1^,*),  which  fit  the  field  data 
better  than  the  previous  solution  set  (see  Figure  7),  and  returned  the  improved 
residuals  J3  =  5.6266  x  10-4  and  J3  =  2.1284  x  10-4.  We  also  found  that 
|u°(0,  a:)  —  w1(0,  i’)|2  =  0.0600. 

Using  a}.- k  we  then  optimized  (19)  for  pjk.  The  resulting  u2(t,x )  did  not 
appear  to  improve  the  fit  to  the  field  data,  as  illustrated  by  the  residuals  = 
5.6266  x  10-4  and  j|  =  2.1284  x  10-4.  In  addition,  we  noted  that 
|p°  —  p1 12  =  1.9451  x  10-6,  which  also  suggested  that  the  latest  iteration  did 
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Figure  7:  Computed  results  u°(t,x )  (-  -)  and  m1)!,  x)  (•  •)  vs.  field  data  ( — )  for 
Example  3. 

not  bring  much  improvement.  The  computed  solutions  u2(t,  x)  are  not  shown 
here  since  the  difference  between  ul(t,x)  and  u2(t,x)  are  not  perceptible  to  the 
eye  at  this  resolution. 

Finally,  we  used  pl-k  to  optimize  (14)  for  a|.  ■  k.  The  resulting  fit  to  the 
field  data  did  improve  slightly  over  the  previous  fit,  returning  the  residuals 
Jf  =  4.6647  x  10-4  and  Jf  =  1.4301  x  10-4.  We  also  found  that  |w2(0,a:)  — 
w3(0,i,)|2  =  0.0043.  Even  though  there  was  an  improvement  in  the  residuals, 
there  was  no  perceptible  improvement  visually  in  the  graph  of  u3(t,  x)  compared 
to  w2(f ,  x). 

The  probability  densities  p°(r,  j)  and  distribution  P°(r,  j)  are  shown  in  Fig¬ 
ures  8  and  9  respectively.  The  densities  and  distributions  for  p1  are  not  shown 
here  because  the  difference  between  p°  and  p1  is  not  perceptible  to  the  eye  at 
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this  resolution.  Note  that  for  fixed  r,  the  probability  distribution  P( 7)  suggests 
a  bimodal  population  density,  as  obtained  in  [BF],  The  probability  densities 
p(r,j)  do  not  contribute  to  this  conclusion,  and  in  fact,  limiting  densities  them¬ 
selves  may  not  exist  for  measures  in  the  set  {P{g)  :  g  6  G}  for  general  infinite 
dimensional  families  G.  (see  [BF];  the  distributions,  not  the  densities,  can  be 
approximated  as  we  have  done  in  these  computations.) 


Figure  8:  Probability  densities  p°jk  for  Example  3. 


Figure  9:  Probability  distribution  P°k  for  Example  3. 
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